{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2ae6082-f8ad-4005-a677-c414f4761a80",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "from linearmodels.panel import PanelOLS\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.preprocessing import LabelEncoder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b0db958-ecab-4d65-a238-6ec82976dbae",
   "metadata": {},
   "source": [
    "# 1. Analysis 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b70316f9-adbf-488c-865b-a58f6b0046e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel(r\"\\df_individual.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e5eb4b5-06b9-4ffd-a677-2579f11051de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop rows with missing values for simplicity\n",
    "df = df.dropna(subset=['Loneliness', 'EVI', 'NDVI', 'sex', 'marital', 'age', 'education', 'income', 'homeownership', 'religion', 'year'])\n",
    "\n",
    "# Ensure proper encoding of categorical variables\n",
    "df['sido_1'] = LabelEncoder().fit_transform(df['sido'])\n",
    "df['si_gun_gu_1'] = LabelEncoder().fit_transform(df['si_gun_gu'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53a329c7-de6b-4826-87c6-26cc300aa1c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Verify the grouping variable\n",
    "assert df['si_gun_gu_1'].isnull().sum() == 0, \"Grouping variable has missing values.\"\n",
    "assert len(df['si_gun_gu_1'].unique()) > 1, \"Grouping variable should have more than one unique value.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "449b4123-1929-4486-84df-90b5c3195b79",
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# Function to run mixed effects model\n",
    "def run_mixed_model(data, formula, re_formula):\n",
    "    model = smf.mixedlm(formula, data, groups=data[\"si_gun_gu_1\"], re_formula=re_formula)\n",
    "    result = model.fit()\n",
    "    return result\n",
    "\n",
    "# Mixed models\n",
    "formulas = [\n",
    "    'Loneliness ~ EVI + sex + marital + age + education + income + homeownership + religion + C(year)',\n",
    "    'Loneliness ~ NDVI + sex + marital + age + education + income + homeownership + religion + C(year)'\n",
    "]\n",
    "re_formulas = ['~EVI', '~NDVI']\n",
    "\n",
    "# Run mixed models\n",
    "results_loneliness = [run_mixed_model(df, formulas[i], re_formulas[i]) for i in range(2)]\n",
    "\n",
    "# Display the results\n",
    "for i, result in enumerate(results_loneliness):\n",
    "    print(f\"Model {i+1}\")\n",
    "    print(result.summary())\n",
    "    print(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fc87636-0056-4af7-bea1-e306054106d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to extract coefficients and CIs\n",
    "def extract_coefs(results, var):\n",
    "    coefs = results.params[var]\n",
    "    ci = results.conf_int().loc[var]\n",
    "    ci_lower = ci[0]\n",
    "    ci_upper = ci[1]\n",
    "    return coefs, ci_lower, ci_upper\n",
    "\n",
    "# Extract for each model\n",
    "variables = ['EVI', 'NDVI']\n",
    "data_for_plot = []\n",
    "for i, var in enumerate(variables):\n",
    "    coefs, ci_lower, ci_upper = extract_coefs(results_loneliness[i], var)\n",
    "    data_for_plot.append([var, f'Model {i+1}', coefs, ci_lower, ci_upper])\n",
    "\n",
    "df_plot = pd.DataFrame(data_for_plot, columns=['Variable', 'Model', 'Coef', 'CI Lower', 'CI Upper'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc5a9155-e472-4846-b1fa-975ce10c65b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "plt.figure(figsize=(12, 8))\n",
    "sns.set(style=\"whitegrid\")\n",
    "\n",
    "# Variables for adjusting positions\n",
    "offset = 0.2\n",
    "positions = {'EVI': -offset, 'NDVI': 0, }\n",
    "colors = {'EVI': 'blue', 'NDVI': 'green', }\n",
    "\n",
    "# Plotting each variable\n",
    "for var in variables:\n",
    "    subset = df_plot[df_plot['Variable'] == var]\n",
    "    x_positions = range(len(subset))\n",
    "    adjusted_positions = [x + positions[var] for x in x_positions]\n",
    "    \n",
    "    plt.errorbar(adjusted_positions, subset['Coef'], yerr=[subset['Coef'] - subset['CI Lower'], subset['CI Upper'] - subset['Coef']],\n",
    "                 fmt='o', label=var, capsize=5, markersize=18, color=colors[var])\n",
    "    \n",
    "    # Annotate coefficient values slightly above the marker\n",
    "    for j in range(len(subset)):\n",
    "        plt.text(adjusted_positions[j], subset['Coef'].iloc[j] + 0.08, f'{subset[\"Coef\"].iloc[j]:.3f}', \n",
    "                 ha='center', va='bottom', fontsize=24, color=colors[var])\n",
    "\n",
    "# Adding the dashed line at 0\n",
    "plt.axhline(0, color='red', linewidth=2, linestyle='--')\n",
    "\n",
    "# Customizing the plot\n",
    "custom_labels = ['DV: Loneliness']\n",
    "plt.xticks(range(len(custom_labels)), custom_labels, size=22)\n",
    "plt.yticks(size=16)\n",
    "\n",
    "# Increasing legend text size\n",
    "plt.legend(title='Variable', fontsize='xx-large', title_fontsize='xx-large')\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ee2d8f7-a8d4-457f-ae2a-564bcd195864",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "371dcb69-cc84-4c02-9a8c-237a7bcb5d72",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mixed models for Distrust, Exploitation, and Indifference\n",
    "additional_formulas = [\n",
    "    'Distrust ~ Loneliness + sex + marital + age + education + income + homeownership + religion + C(year)',\n",
    "    'Exploitation ~ Loneliness + sex + marital + age + education + income + homeownership + religion + C(year)',\n",
    "    'Indifference ~ Loneliness + sex + marital + age + education + income + homeownership + religion + C(year)'\n",
    "]\n",
    "\n",
    "results_additional = [run_mixed_model(df, formula, '~NDVI') for formula in additional_formulas]\n",
    "\n",
    "# Display the results\n",
    "for i, result in enumerate(results_additional):\n",
    "    print(f\"Model {i+1}\")\n",
    "    print(result.summary())\n",
    "    print(\"\\n\")\n",
    "\n",
    "# Extract and plot coefficients for these models as done before.\n",
    "data_for_additional_plot = []\n",
    "for i, var in enumerate(['Loneliness']*3):\n",
    "    coefs, ci_lower, ci_upper = extract_coefs(results_additional[i], var)\n",
    "    data_for_additional_plot.append([var, f'Model {i+1}', coefs, ci_lower, ci_upper])\n",
    "\n",
    "df_additional_plot = pd.DataFrame(data_for_additional_plot, columns=['Variable', 'Model', 'Coef', 'CI Lower', 'CI Upper'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa7afe67-49f6-4111-b08b-718503dfdef4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "plt.figure(figsize=(12, 8))\n",
    "sns.set(style=\"whitegrid\")\n",
    "\n",
    "# Plotting each variable\n",
    "for var in ['Loneliness']:\n",
    "    subset = df_additional_plot[df_additional_plot['Variable'] == var]\n",
    "    x_positions = range(len(subset))\n",
    "    \n",
    "    plt.errorbar(x_positions, subset['Coef'], yerr=[subset['Coef'] - subset['CI Lower'], subset['CI Upper'] - subset['Coef']],\n",
    "                 fmt='o', label=var, capsize=5, markersize=18, color='blue')\n",
    "    \n",
    "    # Annotate coefficient values slightly above the marker\n",
    "    for j in range(len(subset)):\n",
    "        plt.text(x_positions[j], subset['Coef'].iloc[j] + 0.005, f'{subset[\"Coef\"].iloc[j]:.3f}', \n",
    "                 ha='center', va='bottom', fontsize=24, color='blue')\n",
    "\n",
    "# Adding the dashed line at 0\n",
    "plt.axhline(0, color='red', linewidth=2, linestyle='--')\n",
    "\n",
    "# Customizing the plot\n",
    "custom_labels = ['Loneliness', 'Loneliness', 'Loneliness']\n",
    "plt.xticks(range(len(custom_labels)), custom_labels, size=22)\n",
    "plt.yticks(size=16)\n",
    "\n",
    "# Increasing legend text size\n",
    "#plt.legend(title='Variable', fontsize='large', title_fontsize='large')\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5b1369b-1b96-4e46-9d36-fa64ba7624d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the plot size and style\n",
    "fig, axes = plt.subplots(1, 2, figsize=(24, 10))\n",
    "sns.set(style=\"whitegrid\")\n",
    "\n",
    "# Variables for adjusting positions\n",
    "offset = 0.3\n",
    "variables = ['EVI', 'NDVI']  # Define the variables to be plotted\n",
    "positions = {'EVI': -offset, 'NDVI': +offset}\n",
    "colors = {'EVI': 'red', 'NDVI': 'blue'}\n",
    "\n",
    "# Adjust the offset for annotations to avoid overlap\n",
    "annotation_offsets = {'EVI': 0.2, 'NDVI': -0.2}\n",
    "\n",
    "# Plot 1: Loneliness with EVI, NDVI (x and y switched)\n",
    "for var in variables:\n",
    "    subset = df_plot[df_plot['Variable'] == var]\n",
    "    x_positions = range(len(subset))\n",
    "    adjusted_positions = [x + positions[var] for x in x_positions]\n",
    "    \n",
    "    # Switch x_positions and subset['Coef'] in errorbar\n",
    "    axes[0].errorbar(subset['Coef'], adjusted_positions, \n",
    "                     xerr=[subset['Coef'] - subset['CI Lower'], subset['CI Upper'] - subset['Coef']],\n",
    "                     fmt='o', label=var, capsize=11, markersize=22, color=colors[var])\n",
    "    \n",
    "    # Annotate coefficient values (adjusted for switched axes)\n",
    "    for j in range(len(subset)):\n",
    "        axes[0].text(subset['Coef'].iloc[j] + 10 + annotation_offsets[var], \n",
    "                     adjusted_positions[j], \n",
    "                     f'{subset[\"Coef\"].iloc[j]:.3f}', ha='center', va='bottom', \n",
    "                     fontsize=32, color=colors[var])\n",
    "\n",
    "# Adding the dashed line at 0 (vertical now)\n",
    "axes[0].axvline(0, color='red', linewidth=2, linestyle='--')\n",
    "\n",
    "# Customizing the plot\n",
    "axes[0].set_ylim(-1, 1)  # Adjust y-axis limits\n",
    "axes[0].set_yticks([0])  # Position tick in the center\n",
    "axes[0].set_yticklabels(['DV: Loneliness'], size=32, ha='center')\n",
    "axes[0].tick_params(axis='x', labelsize=30)  # Adjust for x-axis now\n",
    "\n",
    "# Increasing legend text size\n",
    "axes[0].legend(title='Variable', fontsize='xx-large', title_fontsize='xx-large')\n",
    "axes[0].grid(True)\n",
    "\n",
    "# Title for left plot\n",
    "axes[0].set_title('(a) Vegetation Index and Loneliness', fontsize=39, pad=30)\n",
    "\n",
    "# Plot 2: Loneliness with NDVI for Distrust, Exploitation, and Indifference (x and y switched)\n",
    "for var in ['Loneliness']:\n",
    "    subset = df_additional_plot[df_additional_plot['Variable'] == var]\n",
    "    x_positions = range(len(subset))\n",
    "    \n",
    "    # Switch x_positions and subset['Coef'] in errorbar\n",
    "    axes[1].errorbar(subset['Coef'], x_positions, \n",
    "                     xerr=[subset['Coef'] - subset['CI Lower'], subset['CI Upper'] - subset['Coef']],\n",
    "                     fmt='o', label=var, capsize=11, markersize=22, color='green')\n",
    "    \n",
    "    # Annotate coefficient values (adjusted for switched axes)\n",
    "    for j in range(len(subset)):\n",
    "        axes[1].text(subset['Coef'].iloc[j] + 0.015, \n",
    "                     x_positions[j], \n",
    "                     f'{subset[\"Coef\"].iloc[j]:.3f}', ha='right', va='bottom', \n",
    "                     fontsize=32, color='green')\n",
    "\n",
    "# Adding the dashed line at 0 (vertical now)\n",
    "axes[1].axvline(0, color='red', linewidth=2, linestyle='--')\n",
    "\n",
    "# Customizing the plot\n",
    "custom_labels = ['Model 1\\n(DV: Distrust)', 'Model 2\\n(DV: Exploitation)', 'Model 3\\n(DV: Indifference)']\n",
    "axes[1].set_yticks(range(len(custom_labels)))\n",
    "axes[1].set_yticklabels(custom_labels, size=28)\n",
    "axes[1].tick_params(axis='x', labelsize=30)\n",
    "\n",
    "# Increasing legend text size and moving it to the bottom right\n",
    "axes[1].legend(title='Variable', fontsize='xx-large', title_fontsize='xx-large', loc='upper left', bbox_to_anchor=(0, 1))\n",
    "axes[1].grid(True)\n",
    "\n",
    "# Title for right plot\n",
    "axes[1].set_title('(b) Loneliness and Fragmentation', fontsize=39, pad=30)\n",
    "\n",
    "plt.tight_layout(pad=5.0)  # Increase space between plots\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2502808-7e65-4308-905d-789018d728fd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b5727a67-e96d-4928-9e44-5fb0d200171f",
   "metadata": {},
   "source": [
    "# 2. Analysis 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9f58616-99dc-4cd1-9aae-8727920042f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel(r\"C:\\df_district.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0248a28b-79c2-450a-af6a-16e8dc8dda0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Encode si_gun_gu (similar to 'encode si_gun_gu' in Stata)\n",
    "df['region_1'], _ = pd.factorize(df['si_gun_gu'])\n",
    "\n",
    "# Step 2: Save the 'year' column separately before creating dummies\n",
    "df['year_original'] = df['year']\n",
    "\n",
    "# Step 3: Create year dummies (similar to 'i.year' in Stata)\n",
    "df = pd.get_dummies(df, columns=['year'], prefix='year', drop_first=True)\n",
    "\n",
    "# Step 4: Set the MultiIndex for panel data (region_1 for entity and year for time)\n",
    "df = df.set_index(['region_1', 'year_original'])\n",
    "\n",
    "# Function to run fixed effects regression\n",
    "def run_fe_regression(data, dependent_var, independent_vars):\n",
    "    y = data[dependent_var]\n",
    "    X = data[independent_vars]\n",
    "    \n",
    "    # Add a constant (intercept)\n",
    "    X = sm.add_constant(X)\n",
    "    \n",
    "    # Run the fixed-effects model\n",
    "    model = PanelOLS(y, X, entity_effects=True)  # entity_effects=True for fixed effects\n",
    "    results = model.fit(cov_type='clustered', cluster_entity=True)\n",
    "    \n",
    "    return results\n",
    "\n",
    "# Define independent variables for the models\n",
    "independent_vars_base = ['inequality','income', 'homeownership','education', 'age', 'sex', 'urbanization'] + [col for col in df.columns if col.startswith('year_')]\n",
    "\n",
    "# Run three models for EVI, NDVI, and LAI\n",
    "results_evi = run_fe_regression(df, 'turn_out', ['EVI'] + independent_vars_base)\n",
    "results_ndvi = run_fe_regression(df, 'turn_out', ['NDVI'] + independent_vars_base)\n",
    "\n",
    "# Step 5: Extract coefficients for EVI, NDVI, and LAI\n",
    "coeffs = {\n",
    "    'EVI': results_evi.params['EVI'],\n",
    "    'NDVI': results_ndvi.params['NDVI'],\n",
    "}\n",
    "\n",
    "errors = {\n",
    "    'EVI': results_evi.std_errors['EVI'],\n",
    "    'NDVI': results_ndvi.std_errors['NDVI'],\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "569fdd77-ab4a-4972-ad5e-cecd8f0e6644",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 6: Plot coefficients with additional formatting and adjusted size\n",
    "fig, ax = plt.subplots(figsize=(12, 8))  # Set plot size to 12x8 inches\n",
    "\n",
    "# Define colors and markers for each variable\n",
    "colors = ['red', 'blue']\n",
    "markers = ['o', 's']  # Different shapes: circle, square\n",
    "\n",
    "# Define the x positions (reduce the space between the plots)\n",
    "x_positions = [12, 10]  # Move NDVI closer to EVI\n",
    "\n",
    "# Multipliers for 95% and 90% confidence intervals\n",
    "ci_multiplier_95 = 1.96\n",
    "ci_multiplier_90 = 1.645\n",
    "\n",
    "# Errors for 95% and 90% confidence intervals\n",
    "errors_95ci = {key: ci_multiplier_95 * errors[key] for key in errors.keys()}\n",
    "errors_90ci = {key: ci_multiplier_90 * errors[key] for key in errors.keys()}\n",
    "\n",
    "# Plot each variable with 95% and 90% CIs\n",
    "for i, (key, color, marker) in enumerate(zip(coeffs.keys(), colors, markers)):\n",
    "    # Plot 90% CI (narrower confidence interval)\n",
    "    ax.errorbar(x_positions[i], coeffs[key], yerr=errors_90ci[key], fmt=marker, color=color, capsize=3, \n",
    "                markersize=17, label=f'{key} (90% CI)', alpha=0.6, linewidth=4)\n",
    "    \n",
    "    # Plot 95% CI (wider confidence interval)\n",
    "    ax.errorbar(x_positions[i], coeffs[key], yerr=errors_95ci[key], fmt=marker, color=color, capsize=5, \n",
    "                markersize=17, label=f'{key} (95% CI)', alpha=1, linewidth=1)\n",
    "\n",
    "    # Add coefficient values as text\n",
    "    ax.text(x_positions[i], coeffs[key] + (0.01 if coeffs[key] > 0 else -0.01), f'{coeffs[key]:.3f}', \n",
    "            ha='center', color=color, fontsize=18)\n",
    "\n",
    "# Set title and labels\n",
    "ax.set_title('DV: Turnout  ', fontsize=18, weight='bold')  # Title with gaps and bold font\n",
    "ax.set_ylabel('', fontsize=12)\n",
    "\n",
    "# Increase xtick and ytick size\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "\n",
    "# Set the x-ticks to display variable names\n",
    "ax.set_xticks(x_positions)\n",
    "ax.set_xticklabels(['EVI', 'NDVI'], fontsize=20)\n",
    "\n",
    "# Horizontal line at y=0 for reference\n",
    "ax.axhline(0, color='red', linestyle='--')\n",
    "plt.grid(True)\n",
    "\n",
    "# Show plot without legend\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d334b27c-9298-4324-ad3c-0e5d4f237df2",
   "metadata": {},
   "outputs": [],
   "source": [
    "evi_min, evi_max = df['EVI'].min(), df['EVI'].max()\n",
    "ndvi_min, ndvi_max = df['NDVI'].min(), df['NDVI'].max()\n",
    "\n",
    "# Generate predicted values based on min and max of EVI and NDVI\n",
    "evi_range = np.linspace(evi_min, evi_max, 100)\n",
    "ndvi_range = np.linspace(ndvi_min, ndvi_max, 100)\n",
    "\n",
    "# Assume the model is simple linear, e.g., pred = intercept + coeff * variable\n",
    "#intercept = 0.2  # Replace with your actual intercept\n",
    "\n",
    "# Predicted values\n",
    "pred_evi = intercept + coeffs['EVI'] * evi_range\n",
    "pred_ndvi = intercept + coeffs['NDVI'] * ndvi_range\n",
    "\n",
    "# Create subplots (1 row, 2 columns)\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))\n",
    "\n",
    "# Subplot 1: Coefficient Plot\n",
    "colors = ['red', 'blue']\n",
    "markers = ['o', 'o']\n",
    "\n",
    "# Bring the EVI and NDVI plots closer together\n",
    "x_positions = [1, 1.00001]  # Make the positions for EVI and NDVI closer\n",
    "\n",
    "# Multipliers for 95% and 90% confidence intervals\n",
    "ci_multiplier_95 = 1.96\n",
    "ci_multiplier_90 = 1.645\n",
    "\n",
    "# Errors for 95% and 90% confidence intervals\n",
    "errors_95ci = {key: ci_multiplier_95 * errors[key] for key in errors.keys()}\n",
    "errors_90ci = {key: ci_multiplier_90 * errors[key] for key in errors.keys()}\n",
    "\n",
    "# Plot each variable with 95% and 90% CIs on ax1\n",
    "for i, (key, color, marker) in enumerate(zip(coeffs.keys(), colors, markers)):\n",
    "    \n",
    "    # Plot 95% CI (wider confidence interval)\n",
    "    ax1.errorbar(x_positions[i], coeffs[key], yerr=errors_95ci[key], fmt=marker, color=color, capsize=5, \n",
    "                markersize=17, label=f'{key} (95% CI)', alpha=1, linewidth=1)\n",
    "\n",
    "    # Add coefficient values as text\n",
    "    ax1.text(x_positions[i], coeffs[key] + (0.02 if coeffs[key] > 0 else -0.01), f'{coeffs[key]:.3f}', \n",
    "            ha='center', color=color, fontsize=20)\n",
    "\n",
    "# Set title and labels for ax1\n",
    "ax1.set_title('(a) Vegetation Index and Turnout  ', fontsize=24, pad=32)  # Increase padding\n",
    "ax1.set_ylabel('', fontsize=12)\n",
    "ax1.tick_params(axis='x', labelsize=18)\n",
    "ax1.tick_params(axis='y', labelsize=14)\n",
    "ax1.set_xticks(x_positions)\n",
    "ax1.set_xticklabels(['EVI', 'NDVI'], fontsize=20)\n",
    "ax1.axhline(0, color='red', linestyle='--')\n",
    "ax1.grid(True)  # Add grid for ax1\n",
    "\n",
    "# Subplot 2: Predicted Value Plot\n",
    "ax2.plot(evi_range, pred_evi, label='Predicted Turnout (EVI)', color='red', linewidth=2)\n",
    "ax2.plot(ndvi_range, pred_ndvi, label='Predicted Turnout (NDVI)', color='blue', linewidth=2, linestyle='--')  # Dotted line for NDVI\n",
    "\n",
    "# Annotate min and max values for EVI\n",
    "ax2.text(evi_min, pred_evi[0], f'{pred_evi[0]:.2f}', ha='right', fontsize=20, color='red')\n",
    "ax2.text(evi_max, pred_evi[-1], f'{pred_evi[-1]:.2f}', ha='left', fontsize=20, color='red')\n",
    "\n",
    "# Annotate min and max values for NDVI\n",
    "ax2.text(ndvi_min, pred_ndvi[0], f'{pred_ndvi[0]:.2f}', ha='right', fontsize=20, color='blue')\n",
    "ax2.text(ndvi_max, pred_ndvi[-1], f'{pred_ndvi[-1]:.2f}', ha='left', fontsize=20, color='blue')\n",
    "\n",
    "# Set labels and title for ax2\n",
    "ax2.set_title('(b) Predicted Turnout from Min to Max of\\n EVI and NDVI', fontsize=24, pad=20)\n",
    "ax2.set_xlabel('', fontsize=14)\n",
    "ax2.set_ylabel('', fontsize=14)\n",
    "ax2.tick_params(axis='x', labelsize=18)  # Increase xtick size\n",
    "ax2.tick_params(axis='y', labelsize=14)  # Increase ytick size\n",
    "ax2.grid(True)  # Add grid for ax2\n",
    "\n",
    "# Add legends to both plots\n",
    "#ax1.legend()\n",
    "ax2.legend(fontsize=16)  # Increase legend size for Plot 2\n",
    "\n",
    "# Show plot\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a5b4797-dcc3-4c85-b26a-a0e49bb5fc23",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
